[Ramses] Add analytical gravity#1673
Conversation
|
Thanks @nbrucy for opening this PR! You can do multiple things directly here: Once the workflow completes a message will appear displaying informations related to the run. Also the PR gets automatically reviewed by gemini, you can: |
Summary of ChangesHello @nbrucy, I'm Gemini Code Assist1! I'm currently reviewing this pull request and will post my feedback shortly. In the meantime, here's a summary to help you and other reviewers quickly get up to speed! This pull request introduces the capability to simulate analytical gravity within the Ramses model. This feature allows for the inclusion of predefined gravitational potentials, such as a point mass, directly into the hydrodynamic simulations. It provides a flexible way to model systems where the gravitational field is known analytically, simplifying setup for specific astrophysical problems like accretion disks, without requiring a full N-body or Poisson solver. The changes span from core C++ solver logic to a new Python example, ensuring both functionality and demonstrability. Highlights
Changelog
Activity
Using Gemini Code AssistThe full guide for Gemini Code Assist can be found on our documentation page, here are some quick tips. Invoking Gemini You can request assistance from Gemini at any point by creating a comment using either
Customization To customize Gemini Code Assist for GitHub experience, repository maintainers can create a configuration file and/or provide a custom code review style guide (such as PEP-8 for Python) by creating and adding files to a Limitations & Feedback Gemini Code Assist may make mistakes. Please leave feedback on any instances where its feedback is incorrect or counter productive. You can react with 👍 and 👎 on @gemini-code-assist comments. If you're interested in giving your feedback about your experience with Gemini Code Assist for Github and other Google products, sign up here. You can also get AI-powered code generation, chat, as well as code reviews directly in the IDE at no cost with the Gemini Code Assist IDE Extension. Footnotes
|
![gemini-code-assist[bot]](https://avatars.githubusercontent.com/in/956858?s=60&v=4){kind=small}
There was a problem hiding this comment.
Code Review
The pull request introduces analytical gravity to the Ramses model, specifically a point mass analytical gravity. This involves adding new configuration structures, updating the solver graph to include a ComputeAnalyticalGravity node, and integrating the gravitational force into the time integration step. The changes also include a Python example demonstrating the usage of this new feature. The code appears to be well-structured and follows the existing patterns within the shamrock framework.
| if cs_expo == 0.5: | ||
| omega = point_mass_Gmass / (rc_soft * rc_soft * rs_soft) | ||
| omega -= (4.0 - cs_expo) * (cs * cs / rc_soft * rc_soft) |
There was a problem hiding this comment.
The calculation for omega when cs_expo == 0.5 seems to have a potential issue with rc_soft * rc_soft appearing twice in the denominator and numerator, respectively. This could lead to incorrect physical results. It looks like a typo where rc_soft * rc_soft might have been intended as rc_soft**2 or rc_soft in the denominator, and rc_soft**2 in the numerator. Please verify the physical formula for this specific case.
omega = point_mass_Gmass / (rc_soft**2 * rs_soft)
omega -= (4.0 - cs_expo) * (cs**2 / rc_soft**2)
omega = max(omega, 0)| rs_soft = np.sqrt(x**2 + y**2 + z**2 + gravity_softening**2) | ||
| x_soft = x * (rs_soft / rs) | ||
| y_soft = y * (rs_soft / rs) |
There was a problem hiding this comment.
The calculation of x_soft and y_soft involves rs_soft / rs. If rs is zero (e.g., at the exact center of the coordinate system), this will result in a division by zero error. While gravity_softening is used to prevent rs_soft from being zero, rs itself can still be zero. Consider adding a small epsilon to rs or handling this edge case explicitly to prevent potential runtime errors.
rs = np.sqrt(x**2 + y**2 + z**2)
rs_soft = np.sqrt(x**2 + y**2 + z**2 + gravity_softening**2)
# Add a small epsilon to rs to prevent division by zero
rs_safe = rs + 1e-10 # or a more appropriate small value
x_soft = x * (rs_soft / rs_safe)
y_soft = y * (rs_soft / rs_safe)| for (u32 d = 0; d < shambase::VectorProperties<Tvec>::dimension; d++) { | ||
| ekin_old += rhov_old[d] * rhov_old[d]; | ||
| } | ||
| ekin_old *= 0.5 / rho[id_a]; |
There was a problem hiding this comment.
The calculation of ekin_old and ekin_new involves division by rho[id_a]. If rho[id_a] can be zero or very close to zero, this could lead to division by zero errors or numerical instability. Consider adding a check for rho[id_a] being non-zero or adding a small epsilon to the denominator to prevent such issues.
Tscal rho_safe = rho[id_a] + 1e-10; // Add a small epsilon
ekin_old *= 0.5 / rho_safe;| for (u32 d = 0; d < shambase::VectorProperties<Tvec>::dimension; d++) { | ||
| ekin_new += rhov[id_a][d] * rhov[id_a][d]; | ||
| } | ||
| ekin_new *= 0.5 / rho[id_a]; |
There was a problem hiding this comment.
Similar to the previous comment, the calculation of ekin_new involves division by rho[id_a]. This also needs a check for zero or near-zero values of rho[id_a] to ensure numerical stability and prevent division by zero errors.
Tscal rho_safe = rho[id_a] + 1e-10; // Add a small epsilon
ekin_new *= 0.5 / rho_safe;|
|
||
| scale_fact = 1 / (cell_size * base) | ||
| cfg.set_scale_factor(scale_fact) | ||
| cfg.set_eos_gamma(1.4) # Can we have something which is NOT an EOS? |
There was a problem hiding this comment.
The comment Can we have something which is NOT an EOS? suggests uncertainty about the eos_gamma setting. If this parameter is not always related to an Equation of State, it might be beneficial to rename it or add a more descriptive comment explaining its broader applicability within the shamrock framework.
| (0, 0, 0), (cell_size, cell_size, cell_size), (base, base, base) | ||
| ) # What is this doing? |
There was a problem hiding this comment.
The comment What is this doing? indicates a lack of clarity regarding the make_base_grid function. It would be helpful to add a more detailed explanation of its purpose and parameters, especially for future maintainers or users of this example script.
| using Tscal = shambase::VecComponent<Tvec>; | ||
| GravityMode gravity_mode = NoGravity; | ||
| bool analytical_gravity = false; // whether to use an external analytical gravity | ||
| Tscal tol = 1e-6; |
There was a problem hiding this comment.
The analytical_gravity boolean was removed from GravityConfig but its comment // whether to use an external analytical gravity was also removed. It would be good to ensure that the purpose of GravityConfig and AnalyticalGravityConfig are clearly delineated, and that no important contextual information is lost.
| for (u32 d = 0; d < dim; d++) { | ||
| r_squared += r_vec[d] * r_vec[d]; | ||
| } | ||
|
|
There was a problem hiding this comment.
The pow function in C++ operates on double by default. Using -1.5f suggests a float literal, which might lead to an implicit conversion and potential precision loss if Tscal is f64. It's generally safer to use pow(base, exponent) where base and exponent have the same type as Tscal to avoid unexpected behavior or to explicitly cast the exponent to Tscal.
Tscal one_over_r_cube_soft = pow(r_squared + config.epsilon_softening, static_cast<Tscal>(-1.5));
|
Wonderfull ! |
There was a problem hiding this comment.
can you
git mv exemples/ramses/godunov_disc.py examples/ramses/run_godunov_disc.py
to account for the new structure of the folder + enable the test to be run by the CI for the doc
There was a problem hiding this comment.
Ok, but for now it's still in debug mode :)
There was a problem hiding this comment.
i mean it's not like there is a choice since the original folder is gone on the main branch 😅
If you don't want to run the script yet just to the same command but do not rename the file.
| }; | ||
|
|
||
| template<class Tvec> | ||
| struct AnalyticalGravityConfig { |
There was a problem hiding this comment.
Regarding the structure of the config I would opt for something closer to the one in SPH
Or just use src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp directly.
But this can be done for a latter PR. As you prefer
|
|
||
| AnalyticalGravityMode analytical_gravity_mode = POINTMASS; | ||
| // parameters for analytical gravity can be added here | ||
| Tscal point_mass_GM = 1; |
There was a problem hiding this comment.
Since G is set by the unit system i would just provide the mass directly
| Tscal one_over_r_cube_soft = pow( | ||
| r_squared + config.epsilon_softening * config.epsilon_softening, -1.5f); |
There was a problem hiding this comment.
| Tscal one_over_r_cube_soft = pow( | |
| r_squared + config.epsilon_softening * config.epsilon_softening, -1.5f); | |
| Tscal one_over_r_cube_soft = sycl::pow( | |
| r_squared + config.epsilon_softening * config.epsilon_softening, Tscal(1.5)); |
Workflow reportworkflow report corresponding to commit 4910b0a Light CI is enabled. This will only run the basic tests and not the full tests. Pre-commit check reportSome failures were detected in base source checks checks. ❌ trailing-whitespaceSuggested changesDetailed changes :diff --git a/src/shammodels/ramses/src/Solver.cpp b/src/shammodels/ramses/src/Solver.cpp
index dbd90f9c..4e79aed0 100644
--- a/src/shammodels/ramses/src/Solver.cpp
+++ b/src/shammodels/ramses/src/Solver.cpp
@@ -25,8 +25,8 @@
#include "shammodels/ramses/SolverConfig.hpp"
#include "shammodels/ramses/modules/AMRGridRefinementHandler.hpp"
#include "shammodels/ramses/modules/BlockNeighToCellNeigh.hpp"
-#include "shammodels/ramses/modules/ComputeAnalyticalGravity.hpp"
#include "shammodels/ramses/modules/ComputeAMRLevel.hpp"
+#include "shammodels/ramses/modules/ComputeAnalyticalGravity.hpp"
#include "shammodels/ramses/modules/ComputeCFL.hpp"
#include "shammodels/ramses/modules/ComputeCellAABB.hpp"
#include "shammodels/ramses/modules/ComputeCoordinates.hpp"
@@ -989,9 +989,8 @@ void shammodels::basegodunov::Solver<Tvec, TgridVec>::init_solver_graph() {
solver_sequence.push_back(
std::make_shared<decltype(node_analytical_gravity)>(
std::move(node_analytical_gravity)));
-
- }
-
+ }
+
if (solver_config.amr_mode.need_amr_level_compute()) { // compute block amr level in patch
modules::ComputeAMRLevel<TgridVec> node_amr_level{};
node_amr_level.set_edges(
|
2 participants
No description provided.